#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "sparse_matrix.h"


//////////////////////////////// sparse matrix operation /////////////////////
/**
 * Sparse matrix add operation. The sparse matrix is represent by matrix J
 * and S.  Where J is the column index whose size is dim x MAXCOUPLE, and
 * the entries are double type and it's meaning is the entry of represented
 * Matrix. We assume that the vale of MAXCOUPLE is big enough and the first
 * column of J represent the number of nonzeros in current row.  so the
 * first column is initialized to be 1, it says that the pointer is at the
 * first position.
 */
void sparse_add(int **J0, double **S0, int row, int col, double val)
{
  int pos;
  // find if value col is in the row's row of matrix J
  for(pos=1; pos <= J0[row][0]; pos++){
    if(J0[row][pos] == col) break;
  }
  // whether J exist or not, i is the current position to insert
  // since S is initialized to be zero, add on is safe, also the rewrite of
  // J[row][i] is also safe!! This implementation is wonderful!
  J0[row][pos] = col;
  S0[row][pos] += val;
  J0[row][0] = (pos > J0[row][0])?pos:J0[row][0];
}

/**
 * treat the sparsity pattern from matrix J0 and S0 the row compress
 * storage format is gotten : I,J,S we assume the memory is reserved by
 * calling functions.  for ease of matrix iteration method, the diagnal
 * entries are put at the beginning of each row.
 */
void sparse_compress(int **J0, double **S0, int *I, int *J, double *S, int dim)
{
  int i,j,pos;
  // fill vector I first, whose first element must be 0
  I[0] = 0;
  for(i = 1; i <= dim; i++){
    I[i] = I[i-1] + J0[i-1][0];
  }
  // then run over J0 and S0 to fill J and S
  for(i = 0; i < dim; i++){
    for(j = 1,pos = 1; j <= J0[i][0]; j++){
      // remember that pos=0 is for the diagnal element
      if(J0[i][j] == i){
    J[I[i]] = i;
    S[I[i]] = S0[i][j];
      }
      else{
    J[I[i]+pos] = J0[i][j];
    S[I[i]+pos] = S0[i][j];
    pos++;
      }
    }
  }
}

/**
 * get the specific entry value on position (i,j)
 */
double sparse_entry(int i, int j, int *I, int *J, double *S, int dim)
{
  // for ease, we assume that i and j are not beyond dim
  // but for purpose of later use, we place dim here
  int pos;
  if (i==j){ // the diagnal element is easy to get!
    pos = I[i];
    return S[pos];
  }
  // otherwise the other element
  for(pos = I[i]+1; pos < I[i+1]; pos++)
    if(J[pos] == j) break;
  if(pos == I[i+1])
    return 0.0; // (i,j) are not int the pattern, so just return 0.0
  else
    return S[pos];  
}

/**
 * Solving linear system using various sparse matrix solver.
 * Currently, gauss-sediel works well for elliptic problem.
 */
void GS(int *I, int *J, double *S, double *rhs, double *x, int dim){
  double err = 1.0,res;
  int iter = 1,i,j,l;
  // double *x_old = (double *) calloc(dim,sizeof(double));
  while(err > 1e-12){
    // this is gauss sediel iteration:
    for(i = 0; i < dim; i++){
      res = rhs[i];
      for(l = I[i] + 1; l < I[i+1]; l++)
    res -= x[J[l]]*S[l];
      x[i] = res/S[I[i]];
    }
    //calculate the norm of the residual
    err = 0.0;
    for(i = 0; i < dim; i++){
      res = rhs[i];
      for(l = I[i]; l < I[i+1]; l++)
    res -= x[J[l]]*S[l];
      err += res*res;
    }
    err = sqrt(err);

    // if((iter++ % 100) == 0)
    //  printf("=== Iteration %d: the error is %le ===\n",iter-1, err);
  }
  printf("The G-S iterative solver finished at step %d with residual = %le.\n",iter-1, err);
}
